#ifndef _IMPROVED_PERLIN_NOISE_H_
#define _IMPROVED_PERLIN_NOISE_H_

/*
    Ken Perlin's "Improved Noise"
    http://mrl.nyu.edu/~perlin/Noise/
*/

#include <iostream>
#include <fstream>

namespace ImprovedPerlinNoise 
{
    static int p[] = 
    { 
        151,160,137,91,90,15,
        131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
        190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
        88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
        77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
        102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
        135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
        5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
        223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
        129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
        251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
        49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
        138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
        151,160,137,91,90,15,
        131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
        190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
        88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
        77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
        102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
        135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
        5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
        223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
        129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
        251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
        49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
        138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
    };

    static float g[][3] = 
    {
         1, 1, 0,
        -1, 1, 0,
         1,-1, 0,
        -1,-1, 0,  
         1, 0, 1,
        -1, 0, 1,
         1, 0,-1,
        -1, 0,-1,  
         0, 1, 1,
         0,-1, 1,
         0, 1,-1,
         0,-1,-1,
         1, 1, 0,
         0,-1, 1,
        -1, 1, 0,
         0,-1,-1,
    };

    static float  fade(float t);
    static float  lerp(float t, float a, float b);
    static float  dot(float a, float b, float c, 
                      float x, float y, float z);

    static float  Noise(float3 p);
    static float  Noise(float x, float y, float z);
    static float  Noise(float x, float y, float z, float w);
    static float3 Noise3f(float3 p);

    static float  BrownianMotion   ( float3 p, 
                                     int    octaves    = 4, 
                                     float  lacunarity = 2.0f, 
                                     float  gain       = 0.5f);

    static float3 BrownianMotionf3 ( float3 p, 
                                     int    octaves    = 4, 
                                     float  lacunarity = 2.0f, 
                                     float  gain       = 0.5f);
    static float  Turbulence       ( float3 p, 
                                     int    octaves    = 4, 
                                     float  lacunarity = 2.0f, 
                                     float  gain       = 0.5f);

    static float  Grad(int hash, float x, float y, float z);
    static float  Grad(int hash, float x, float y, float z, float w); 

    void  dump_perm_2d(char *filename);

    static float Noise(float x, float y, float z) 
    {
      int X = (int)floor(x) & 255;                                    // FIND UNIT CUBE THAT
      int Y = (int)floor(y) & 255;                                    // CONTAINS POINT.
      int Z = (int)floor(z) & 255;
      x -= floor(x);                                                  // FIND RELATIVE X,Y,Z
      y -= floor(y);                                                  // OF POINT IN CUBE.
      z -= floor(z);
      float u = fade(x);                                              // COMPUTE FADE CURVES
      float v = fade(y);                                              // FOR EACH OF X,Y,Z.
      float w = fade(z);
      int A = p[X  ]+Y, AA = p[A]+Z, AB = p[A+1]+Z,                   // HASH COORDINATES OF
          B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z;                   // THE 8 CUBE CORNERS,

      return lerp(w, lerp(v, lerp(u, Grad(p[AA  ], x  , y  , z   ),   // AND ADD
                                     Grad(p[BA  ], x-1, y  , z   )),  // BLENDED
                             lerp(u, Grad(p[AB  ], x  , y-1, z   ),   // RESULTS
                                     Grad(p[BB  ], x-1, y-1, z   ))), // FROM  8
                     lerp(v, lerp(u, Grad(p[AA+1], x  , y  , z-1 ),   // CORNERS
                                     Grad(p[BA+1], x-1, y  , z-1 )),  // OF CUBE
                             lerp(u, Grad(p[AB+1], x  , y-1, z-1 ),
                                     Grad(p[BB+1], x-1, y-1, z-1 ))));
    }

    static float Noise(float3 p)
    {
        return Noise(p.x, p.y, p.z);
    }

    // vector Noise
    static float3 Noise3f(float3 p)
    {
        return f3(Noise(p),
                         Noise(p + f3(32, 78, 7)), 
                         Noise(p + f3(123, 11, 96)));
    }

    // fractal sum
    static float BrownianMotion(float3 p, 
                                int    octaves,  
                                float  lacunarity, 
                                float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;	
        for(int i=0; i<octaves; i++) 
        {
	        noiseSum += Noise(p * freq) * amp;

	        freq *= lacunarity;
	        amp  *= gain;
        }
        return noiseSum;
    }

    static float3 BrownianMotionf3(float3 p, 
                                   int    octaves, 
                                   float  lacunarity, 
                                   float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float3 noiseSum = {0.0f, 0.0f, 0.0f};	
        for(int i=0; i<octaves; i++) 
        {
	        noiseSum += Noise3f(p * freq) * amp;

	        freq *= lacunarity;
	        amp  *= gain;
        }
        return noiseSum;
    }

    static float Turbulence(float3 p, int octaves, float lacunarity, float gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;
        for(int i=0; i<octaves; i++) 
        {
	        noiseSum += abs(Noise(p * freq)) * amp;

	        freq *= lacunarity;
	        amp  *= gain;
        }
        return noiseSum;
    }

    static float fade(float t) 
    { 
        return t * t * t * (t * (t * 6 - 15) + 10); 
    }

    static float lerp(float t, float a, float b) 
    { 
        return a + t * (b - a); 
    }

    static float dot(float a, float b, float c, 
                     float x, float y, float z) 
    { 
        return a*x + b*y + c*z; 
    }

    #if 0
    float Grad(int hash, float x, float y, float z) {
      int h = hash & 15;                      // CONVERT LO 4 BITS OF HASH CODE
      float u = h<8 ? x : y;                 // INTO 12 GRADIENT DIRECTIONS.
      float v = h<4 ? y : h==12||h==14 ? x : z;
      return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
    }
    #else
    // use look-up table
    static float Grad(int hash, float x, float y, float z) 
    {
        int h = hash & 15;
        return dot(g[h][0], g[h][1], g[h][2], x, y, z);
    }
    #endif

    static float Grad(int hash, float x, float y, float z, float w) 
    {
        int h = hash & 31;                 // CONVERT LO 5 BITS OF HASH 
                                           // TO 32 GRAD DIRECTIONS.
        float a=y,b=z,c=w;                 // X,Y,Z
        switch (h >> 3) 
        {                                  // OR, DEPENDING ON HIGH ORDER 2 BITS:
            case 1: a=w;b=x;c=y;break;     // W,X,Y
            case 2: a=z;b=w;c=x;break;     // Z,W,X
            case 3: a=y;b=z;c=w;break;     // Y,Z,W
        }
        return ((h&4)==0 ? -a:a) + ((h&2)==0 ? -b:b) + ((h&1)==0 ? -c:c);
    }

    static float Noise(float x, float y, float z, float w)
    {
        int X = (int)floor(x) & 255;                                 // FIND UNIT HYPERCUBE
        int Y = (int)floor(y) & 255;                                 // THAT CONTAINS POINT.
        int Z = (int)floor(z) & 255;
        int W = (int)floor(w) & 255;

        x -= floor(x);                                              // FIND RELATIVE X,Y,Z,W
        y -= floor(y);                                              // OF POINT IN CUBE.
        z -= floor(z);
        w -= floor(w);
        float a = fade(x);                                          // COMPUTE FADE CURVES
        float b = fade(y);                                          // FOR EACH OF X,Y,Z,W.
        float c = fade(z);
        float d = fade(w);

        int A = p[X  ]+Y, AA = p[A]+Z, AB = p[A+1]+Z;               // HASH COORDINATES OF
        int B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z;               // THE 16 CORNERS OF
        int AAA = p[AA]+W, AAB = p[AA+1]+W;                         // THE HYPERCUBE.
        int ABA = p[AB]+W, ABB = p[AB+1]+W;
        int BAA = p[BA]+W, BAB = p[BA+1]+W;
        int BBA = p[BB]+W, BBB = p[BB+1]+W;
                                                                    // INTERPOLATE DOWN.
        return lerp(d, lerp(c,lerp(b,lerp(a,Grad(p[AAA  ], x  , y  , z  , w  ), 
                                            Grad(p[BAA  ], x-1, y  , z  , w  )),
                                     lerp(a,Grad(p[ABA  ], x  , y-1, z  , w  ), 
                                            Grad(p[BBA  ], x-1, y-1, z  , w  ))),

                              lerp(b,lerp(a,Grad(p[AAB  ], x  , y  , z-1, w  ), 
                                            Grad(p[BAB  ], x-1, y  , z-1, w  )),
                                     lerp(a,Grad(p[ABB  ], x  , y-1, z-1, w  ),
                                            Grad(p[BBB  ], x-1, y-1, z-1, w  )))),

                       lerp(c,lerp(b,lerp(a,Grad(p[AAA+1], x  , y  , z  , w-1), 
                                            Grad(p[BAA+1], x-1, y  , z  , w-1)),
                                     lerp(a,Grad(p[ABA+1], x  , y-1, z  , w-1), 
                                            Grad(p[BBA+1], x-1, y-1, z  , w-1))),

                              lerp(b,lerp(a,Grad(p[AAB+1], x  , y  , z-1, w-1), 
                                            Grad(p[BAB+1], x-1, y  , z-1, w-1)),
                                     lerp(a,Grad(p[ABB+1], x  , y-1, z-1, w-1),
                                            Grad(p[BBB+1], x-1, y-1, z-1, w-1)))));
    }

    static void dump_perm_2d(char *filename)
    {
			std::ofstream file("noise.txt");

      int Z = 0;
      for(int Y=0; Y<256; Y++) 
      {
        for(int X=0; X<256; X++) 
        {
          int  A = p[X  ]+Y, AA = p[A]+Z, AB = p[A+1]+Z;
          int  B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z;
          char c[4];
          c[0] = AA;
          c[1] = BA;
          c[2] = AB;
          c[3] = BB;
          file << c[0] << c[1] << c[2] << c[3];
        }
      }
			file.close();
    }
}

#endif